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Abstract 

A computationally-fast inverse design method for nanophotonic struc- 
tures is presented. The method is based on two complementary convex 
optimization problems which modify the dielectric structure and resonant 
field respectively. The design of one- and two-dimensional nanophotonic 
resonators is demonstrated and is shown to require minimal computational 
resources. 
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1 Introduction 

Numerous numerical methods have been devised to solve Maxwell's equations in 
both time[l] and frequency [2^, domains. We refer to these schemes as direct 
solvers, since they compute the electric and magnetic fields based on current 
sources, charge distributions and surrounding dielectric and/or metallic struc- 
tures. While extremely useful in simulating optical components, using direct 
methods to design optical components, especially in two or three dimensions, 
typically requires an extremely time-consuming brute force search in a large 
parameter space [H [H [H [71 [8 . 

On the other hand, an inverse solver would be much more adept in such 
design and optimization problems [9, 10 . In this work, we define the inverse 
problem as that where the electromagnetic field is known, but the surrounding 
structure is not known. The goal in the inverse problem is, then, to find a 
dielectric structure that will produce that specific electromagnetic field profile. 

We show that one can design nanophotonic resonators by specifying the 
electromagnetic field and its desirable characteristics (such as cavity quality 
(Q) factor and/or mode volume) and then using an inverse solver to find the 
corresponding dielectric structure. We show that the inverse method used is 
not only computationally-fast, but is also able to optimize for multiple device 
characteristics and produce multiple resonances, both of which are very difficult 
using direct methods. 



2 Numerical Setup 

We start from the time-harmonic eigenvalue equation 

V X e-^V xH = H (1) 

where e, co and c are the magnetic field, relative permittivity, resonance 
frequency and speed of light respectively. To solve the problem numerically, H 
and e are discretized in space using the standard Yee cell used in finite difference 
methods[l]. Also, the curl operators, since they are linear, are represented by 
the matrix A. Eq. ([T]) can now be written as 

AY Ax = (2) 

where 

A is the discretized curl operator, 

Y = diag(e~^) is the diagonal matrix representing the dielectric structure, 
X is a vector representing and 



2 



In this form, given F, we can solve the direct problem by computing x using an 
eigenvalue solver [1 2 . However, we note that eq. ([2| is also linear in F, which 
allows us, if X is held constant, to solve the inverse problem by expressing eq. ([T]) 
as 

By = d (3) 

where 



B = A' diag(Ax), 
d = ^x, and 



the variable for which we solve. 



Here, diag(Ax) is the matrix with the values of Ax along the main diagonal and 
zeros elsewhere. 



3 Least-Squares Method in ID 
3.1 Least- Squares 
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Figure 1: Inverse design of a one-dimensional structure using the unmodified 
least-squares method. The target field is a sinusoid within a Gaussian envelope. 
The computed dielectric structure (green area) supports a field (red circles) that 
exactly matches the target field (blue line). The entire design process is also 
extremely fast and takes less than 1 second to complete on a generic desktop 
computer. The periodic singularities in the dielectric structure are non-physical 
and will be addressed later in the article. 



The result of applying this method to a simple one-dimensional problem is 
shown in Fig.jl] A generic least-squares solver [11 was used to find the dielectric 
structure, y (green region), that exactly produces the target field, x (blue line), 
using eq. ([2|. Using a generic desktop computer, the solution was obtained in 
less than a second. Then a finite-difference time-domain (FDTD) solver was 
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used to obtain the actual field (red circles) produced by the structure and to 
verify the accuracy of y. 

As expected, Fig. [l] shows that the target field is reproduced exactly by the 
dielectric structure. However, the resulting structure is full of undesireable sin- 
gularities. The rest of the section focuses on producing a well-behaved dielectric 
structure that still reproduces the target field accurately. 

3.2 Regularized Least- Squares 

The simplest way to produce a well-behaved dielectric structure is to add a 
regularization term to our least-squares problem. 



■ B ' 




d 




y = 


_Vvyo_ 



(4) 



which is equivalent to solving the following optimization problem 

minimize \\By - df ^ r]\\y - yof . (5) 

Here yo represents some initial guess for the dielectric structure, that we want 
the values of y to stay close to. 

This method allows us to trade off accuracy in order to keep e close to 
acceptable values, by increasing r]. We chose to constrain e around a constant 
value of eo = 10 and solved the least-squares system for r] = 10~^, 10~^, and 
10~^. The results, each still obtained in under a second, are shown in Fig. [2] 
and illustrate the trade-off between constraining e and accurately reproducing 
H. 

4 Complementary Optimization in ID 

4.1 Motivation for a Complementary Optimization Strat- 
egy 

The fundamental problem in the previous examples is actually not in the meth- 
ods themselves, but in the improper selection of a target field. In fact, it is very 
difficult to select a suitable target resonant field because not every resonant 
mode even has a corresponding dielectric structure that is able to reproduce it. 
Furthermore, it is nearly impossible to select a multi-dimensional field which 
corresponds to a well-behaved, isotropic and discretely- valued e, as would be 
needed for practical structures. 

For this reason, a successful method must be allowed to either modify the 
target field, or specify it completely, in which case the user would only deter- 
mine certain characteristics (e.g. mode-volume, Q-factor) that the target field 
should have. The former strategy is developed in both one and two dimen- 
sions, while the latter strategy is implemented in Section |6] in order to design 
two-dimensional resonators with discrete values of e. 
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Figure 2: Inverse design of one-dimensional structures using the regularized 
least-squares method. The same target field is used as in Fig. [l] and the com- 
putation time remains below 1 second. As the regularization parameter, r^, is 
increased, e is increasingly constrained to a chosen constant value of 10. At the 
same time, the mismatch between target and actual fields increases markedly. 
This illustrates the apparent trade-off between producing reasonable structures 
and accurately reproducing a fixed target field. 
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4.2 Complementary Optimization 



We start with the same target field as in the previous examples but we now 
formulate a method that allows for it to be modified during the design process. 
The formulation chosen is a complementary optimization routine, where we con- 
tinually alternate between modifying e to better fit the field, and then modifying 
the field to better fit e. Here, we use the term "fit" to mean that either e or H 
is solved so that the residual error from eq. ([T]) is minimized. Additionally, both 
iterations are regularized in order to stably approach a solution. This algorithm 
can be summarized as follows. 



choose xo and 
for z = 1, 2, . . . 

minimize \\Bi_iyi ^ r]i\\yi 

|2 



minimize ||Ay^Axi — (,Xi-i\ 



Vi-iW 



(6) 
(7) 



where Yi = diag(?/i), Bi = A - diag(Axi) and di = ^x^. \\AYiAxi — 
used instead of \\AYiAxi — to avoid the trivial Xi = solution and does 

not affect the overall accuracy since x changes very slowly. 
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Figure 3: Inverse design of a one-dimensional structure using the complementary 
optimization method. The target field in Figs [l] and [2] is used as the initial target 
field. The rates of change for both e and H are controlled by regularization 
parameters r]i = 10~^ and r]2 = 10~^ respectively. The 400 iterations used 
to achieve this result took 60 seconds to compute. This method results in a 
well-behaved e that actually produces a field very similar to the original target 
field. Interestingly, the formation of a "steady-state" periodic structure toward 
the sides of the structure has emerged. 



Fig. |3] shows that the complementary optimization algorithm, after 400 iter- 
ations and with the correct choice of regularization parameters rji and 7^2, results 
in a well-behaved structure that is able to closely reproduce the modified tar- 
get field. Numerically, the least-squares problem must now be solved numerous 
times, which increases the computational time needed to around 60 seconds. 
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4.3 Complementary Optimization with Bounded e 



In order to achieve a more practical, discretely-valued dielectric structure, we 
can impose strict upper- and lower-bounds on e. To this end, we modify our 
algorithm as such, 

choose xo and yo 
for z = 1, 2, . . . 

minimize \\Bi_iyi - 

subject to e'l^ < Vi < e']^ (8) 

minimize HAy^Axi - ^Xi^if ^ r]2\\xi - Xi^if . (9) 

In this algorithm, eq. Q is a convex optimization problem[13 . This allows 
us to impose hard constraints on e, which in turn allows us to remove the 
regularization term present in eq. (|6|. The CVX package[14 is used to solve 
eq. ([8|, with each iteration of the algorithm now requiring roughly 1 second of 
computation time. 
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Figure 4: Inverse design of a one-dimensional structure using the complementary 
optimization method with bounded e. The parameters are identical to those 
used to produce Fig. [S] with the exception that only one regularization term is 
now needed (772 = 10~^). The algorithm was run for 100 iterations, which took 
100 seconds. The structure turns out to be almost completely binary- valued and 
looks like a periodic structure with tapered duty cycle. It produces an actual 
field which very closely matches the final target field. 

A nearly binary- valued dielectric structure is obtained in Fig. [4j which ac- 
curately produces the final target field. This is very useful for the design of 
practical structures, since they usually consist of two or three different ma- 
terials at most. Interestingly, although the directly discreteness of e was not 
enforced (since that would make the problem non-convex), a discrete, binary- 
valued structure has still arisen. 
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Figure 5: Inverse design of an "S" resonator using the complementary opti- 
mization method without bounds on e. The design was initialized by specifying 
an initial dielectric structure (e = 1 everywhere) and a resonant field in the 
shape of an "S" . The final dielectric structure was produced after 50 iterations 
which took 90 seconds to complete in total. The final dielectric structure is 
quite unintuive, and yet reproduces the target field surprisingly well. This ex- 
ample demonstrates the versatility of the complementary optimization method 
in producing designs, from very simple specifications, which otherwise could be 
attained only with considerable difficulty. 



5 Complementary Optimization in 2D 

5.1 "S" Resonator 

We now demonstrate that the complementary optimization method is versatile 
and can be scaled to multiple dimensions. To ensure that e is well-behaved we 
use a point-spread function which does not allow e to change at a certain point 
in space without affecting the values surrounding it. 

In order to show that our method can produce complex designs, we choose 
an S-shaped target field which is non-trivial to reproduce. The optimization 
results, using the complementary optimization method from Section |4.2[ are 
shown in Fig. |5j The resulting dielectric structure is continuous, unbounded 
and contains some singularities (white dots), but the final target and actual 
fields match up well. Also, the computational cost remains quite reasonable; 
the 50 iterations needed required only 5 minutes of computation time. The 
resulting structure is completely unintuitive, and illustrates the kind of new 
capabilities offered by the inverse design strategy. Specifically, that a complex, 
intricate structure can be designed just by specifying the shape and frequency 
of a rather simple electromagnetic mode. 
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5.2 Multi-Mode Inverse Design 



The complementary optimization method can also be extended to produce di- 
electric structures with multiple resonances. To do so, multiple initial target 
fields are specified. The dielectric structure is first modified to simultaneously 
fit all target fields using a mult i- objective least-squares method. Then each 
target field is individually modified to fit the structure; and we continue alter- 
nating between optimizing e and H'^^^ in this way. A benefit of this scheme 
is that only the e optimization increases in size, so the design process remains 
computationally tractable, even for several resonant fields. This algorithm can 
be summarized as follows, 

choose ?/o, Xq'\ Xq \ • • • 
for i = 1, 2, . . . 

minimize ^^y, - y,_^\\^ + Y^\\B\i\y^ - d!f\f (10) 

3 

for j = l,2,... 

minimize WAY^Ax^p - + r]2\\x^P - x^P\f. (11) 

The design of an "X" resonator with two perpendicular, degenerate modes is 
shown in Fig. [6) The added complexity increases the total computation time to 
5 minutes for 40 iterations of the algorithm. 



5.3 Design of Waveguiding Devices 

The multi-mode inverse design method can also be applied to the design of 
waveguiding devices such as multiplexers/demultiplexers, waveguide couplers, 
crossbars and dispersion-tailored waveguides. This can be accomplished by 
treating a waveguiding device as a doubly-degenerate resonator at its oper- 
ational frequencies and then enforcing opposite symmetries (even/odd or co- 
sine/sine) in the degenerate modes. A single-to-dual beam waveguide coupler 
designed based on Eq. (10) and (11) is shown in Fig. [7| and motivates how one 



might design other waveguiding devices such as channel-drop filters. In addi- 
tion to macroscopic waveguiding devices, one should also be able to create novel 
periodic waveguiding structures. For instance, by controlling the frequency of a 
particular waveguiding mode at each of its /c- vectors, one could create a waveg- 
uide with a customized dispersion characteristic, which would be useful in the 
design of slow-light waveguides. 
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Figure 6: Inverse design of a doubly-resonant, degenerate "X" resonator pro- 
duced by the complementary optimization method with unbounded e. As in 
Fig. |5] the initial value of e was 1 everywhere. The two initial target fields used 
are only slightly perturbed and are very similar to the two final actual target 
fields. Computationally, this design took 5 minutes to complete and required 40 
iterations. This example shows that the complementary optimization strategy 
can be extended to produce dielectric structures with multiple resonances. Such 
an "X" resonator is useful for polarization-entangled single-photon sources for 
example [15]. 



6 Complementary Optimization with Bounded 
6 in 2D 

6.1 Numerical Method 

We now use the complementary optimization method to design resonators with 
discrete, binary e in two dimensions. At the same time, the algorithm is mod- 
ified so that we no longer specify an initial target field; instead, only an initial 
dielectric structure and the maximum desired mode volume (i.e. mode area in 
2D) are specified. In addition, the optimization process attempts to maximize 
the Q-factor. Such an algorithm now consists of iteratively solving two convex 
optimization problem, 

choose yo 

for i = 1,2,... 

minimize HAY^-iAx^ — + 7^||Fx^|p 

subject to {AxifYi_i{Axi) < A^^de (12) 



minimize \\Biyi - di\\'^ 

subject to e'l^ < Vi < e']^. (13) 
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Figure 7: Inverse design of a single-to-dual beam waveguide coupler using the 
complementary optimization method without bounds on e. Two degenerate 
modes with opposite symmetry (sine and cosine) are used as target fields. Only 
4 iterations (14 seconds) are needed to achieve this solution. This is a simple 
demonstration showing that the complementary optimization method can also 
be extended to design waveguiding devices. 
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In this algorithm, the field optimization, Eq. ( |T2| ), differs significantly from 
Eq. ([9]). First, the eigenvalue equation term ||74y^_iAx^— ^x^lp no longer utilizes 
Xi-i. Also, a new "Fourier-minimizing" term has been added. Here, the 

row vectors of F consist of the field Fourier components that we do not want 
incorporated in the final solution. The motivation for adding this term is to 
design high-Q resonators using planar structures, where the Q-factor is limited 
by out-of-plane losses, and can thus be improved by eliminating field Fourier 
components that are not localized by total internal refiection (components inside 
the light cone) [10 . Lastly, the {Axi)^Yi-i{Axi) term allows us to specify the 
mode area (mode volume in three dimensions), that we desire for our resonant 
field. This works because the two minimization objectives in Eq. (12) will 
generally cause the mode area to be as large as possible, which means we always 
end up with {AiXi)^Yi_i{AiXi) = A^ode- 

The addition of the two new terms in Eq. ( 12 ) signifies that the field iteration 
in our algorithm attempts to do more than to just satisfy Maxwell's equations. 
Rather than only trying to fit the dielectric structure, the resonant field now 
also minimizes some of its Fourier components while working with a limited 
mode area. 

In contrast with the field optimization given by Eq. (12), the structure 
optimization given by Eq. (13) remains identical to the equation in the one- 
dimensional case, Eq. ([8|, and contains no terms related to notions of Fourier 
components or mode area. In fact, the only objective in the e optimization is 
to better fit the field generated from the field optimization. This means that in 
this scheme, the field iteration is leading the structure iteration. In other words, 
Eq. (12) "looks ahead" and gradually adapts itself to become a more desirable 
field, while Eq. (13) just "follows along". For this reason, this algorithm is 
unique from the other complementary optimization algorithms previously pre- 
sented in the article, since in those algorithms both iterations follow one another 
and eventually "meet in the middle" . 

Finally, in this scheme, we chose to limit the degrees of freedom of both 
Xi and yi. Some components of Xi were fixed at non-zero values and were 
not allowed to be modified simply to avoid the Xi = situation. To enhance 
the aesthetics of the resulting dielectric structure, some components of were 
frozen as well, which is useful in cases where one would only like to modify the 
dielectric structure within a waveguide only, for example. 



6.2 Circular Grating Resonator 

Fig. [8] shows the design of a circular grating resonator using the method from 
Section |6.1| The dielectric structure emerged from a very simple choice of 
initial structure, namely a constant e = 12.25 everywhere. The range of e 
was limited to be from 1 to 12.25, the mode area was set to 5 and r] = 10~^. 
Central components of Xi were held constant to ensure that an electric dipole 
in the y-direction was produced in the center of the structure. Additionally, 
the components of e outside a specified circle were held at a constant e = 12.25 
for the duration of the design process. The entire algorithm was run for 40 
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Figure 8: Inverse design of a two-dimensional resonator using the complemen- 
tary optimization method with strict bounds on e. The initial specification is 
very simple and consists an initial dielectric structure (e = 12.25 everywhere), 
the frequency and mode-volume of the resonance field as well as a weighting 
factor, 7^, to avoid leaky field Fourier components. Additionally, the values of 
e are only allowed to be modified within a central circular region and must be 
kept between 1 and 12.25. After 40 iterations which took 7 minutes to complete, 
a discrete structure emerged with excellent match between the predicted {x^o) 
and actual fields. The structure resembles a circular grating with a bowtie-like 
central structure for focusing the resonant energy to a single point. 



iterations and took 7 minutes. 

After 40 iterations, we see that a strong binary-valued structure has formed. 
Interestingly, the central bowtie-like structure has emerged from previous ge- 
neetic optimization methods as well[7j. Note also the extremely small amount 
of information needed to produce this result: only the frequency and mode 
area of the resonant field desired, and a trivial intial e. The ability to produce 
a rather advanced design from such a simple problem specification highlights 
the potential usefulness of inverse methods in the design of novel nanophotonic 
devices. 



6.3 Beam Resonator 

The same approach was used to design a beam resonator as shown in Fig. [9| 
The parameters used in this design are identical to those used for the circular 
resonator, with the exception that the initial dielectric structure now consists of 
an unbroken waveguide, and the region where e can be modified is now confined 
to the center of the waveguide. 

These two, simple modifications produce a vastly different result as is plainly 
seen from Fig. |9] Some characteristics remain, such as the two closely spaced 
holes in the center of the structures which focus the electromagnetic energy 
to a central point. Other features are unique to the beam resonator, such as 
a gradual outward tapering of hole diameters and size ending in a periodic 
"steady-state" structure, in similar fashion to Fig. [3) This tapering is a direct 



result of the Fourier term in Eq. ( 12 ) as it leads to a smooth field profile variation 



and minimization of the components in the light cone that lead to out-of-plane 
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Figure 9: Inverse design of a beam resonator in two dimensions using the com- 
plementary optimization method with bounded e. The initial conditions are 
identical to those for Fig. [Sj except that the initial dielectric structure is an 
unbroken waveguide, and e can only be modified within that waveguide. The 
structure emerged after 40 iterations, which took 5 minutes of computation. 
The bowtie-like structure has reappeared in the center. Interestingly, the effect 



of the Fourier term in Eq. (12), is seen in the outward tapering of the holes. 



losses [TOl [16]. 

7 Conclusions and Outlook 

In summary, we have introduced and demonstrated a complementary optimiza- 
tion method for the inverse design of nanophotonic resonators as well as waveg- 
uiding structures. We have demonstrated that, in two dimensions, this method 
can be used to design non-trivial as well as multi-resonant structures. Fur- 
thermore, we have shown that by constraining the range of available dielectric 
constants, the same underlying method can produce binary- valued, discrete di- 
electric structures. 

The methods presented here can readily be extended into three dimensions. 
In essence, the only hurdle is the size of the convex optimization problem for 
the field optimization. Since a full three-dimensional grid often consists of tens 



of millions of variables, solving the field optimization problem, Eq. (12), using 
methods such as LU factorization is no longer computationally tractable. In- 
stead, iterative methods such as truncated Newton methods must be employed. 
As it turns out, such iterative methods are especially well suited to this comple- 
mentary optimization scheme, since they can take advantage of the information 
in Xi-i in order to solve Xi very quickly. This is very applicable to our algorithm 
since the field changes very little from one iteration to the next. 

Although solving the field optimization problem is numerically challenging 



in three dimensions, solving the structure optimization problem, Eq. (13), in 
three dimensions is much simpler. This is because only structural variations 
in-plane need to be accounted for planar nanophotonic devices. The number 



of variables in Eq. (13) for both two- and three-dimensional problems is then 



14 



effectively the same, and the same numerical techniques can be employed to 
solve both. 

A full three-dimensional inverse design method based on complementary 
optimization would enable the design of novel nanophotonic resonators. For 
instance, a resonator for efficient wavelength-conversion could be designed by 
creating a structure with multiple well-placed resonant frequencies [17 . Further- 
more, effects such as mode overlap and the use of frequency- dependent refractive 
index materials could be taken into account with such a method. 
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